ChannelInitiation Subroutine

public subroutine ChannelInitiation(method, threshold, mask, flowAcc, flowDir, dem, channel)

Define channel cells. Two methods are possible:

  1. area: channel begins where drainage basin exceedes a given area
  2. ask: channel begins where ASk exceedes a given value, where A is the basin area (m2), S is the local slope (-), and k is the basin fractal dimension (default = 1.7)

Arguments

Type IntentOptional Attributes Name
character(len=*), intent(in) :: method
real(kind=float), intent(in) :: threshold
type(grid_integer), intent(in) :: mask
type(grid_integer), intent(in) :: flowAcc
type(grid_integer), intent(in) :: flowDir
type(grid_real), intent(in) :: dem
type(grid_integer), intent(inout) :: channel

Variables

Type Visibility Attributes Name Initial
real, public :: area

( m2 )
real, public :: ask

local value of ASk
integer, public :: i
integer, public :: id
integer, public :: iu
integer, public :: j
integer, public :: jd
integer, public :: ju
real, public :: scale = 1.7
type(grid_real), public :: slope

Source Code

SUBROUTINE ChannelInitiation &
!
(method, threshold, mask, flowAcc, flowDir, dem, channel )


IMPLICIT NONE

!arguments with intent in:
CHARACTER (LEN = *), INTENT(IN) :: method
REAL (KIND = float), INTENT(IN) :: threshold
TYPE (grid_integer), INTENT(IN) :: mask
TYPE (grid_integer), INTENT(IN) :: flowAcc
TYPE (grid_integer), INTENT(IN) :: flowDir
TYPE (grid_real), INTENT(IN) :: dem

!arguments with intent(inout):
TYPE (grid_integer), INTENT(INOUT) :: channel

!local declarations:
INTEGER :: i,j,iu,ju,id,jd
REAL    :: area !! ( m<sup>2</sup> )
REAL    :: scale = 1.7
REAL    :: ask !!local value of  AS<sup>k</sup>
TYPE (grid_real) :: slope

!----------------------end of declarations-------------------------------------

!initialize channel to zero
 DO i = 1, dem % idim
    DO j = 1, dem % jdim
        IF (mask % mat (i,j) /= mask % nodata ) THEN
            channel % mat (i,j) = 0
        END IF
    END DO
 END DO

SELECT CASE (method)
  CASE ('area')
    DO i = 1, mask % idim
      DO j = 1, mask % jdim
        IF (mask % mat (i,j) /= mask % nodata ) THEN
          area = flowAcc % mat (i,j) * CellArea (flowAcc, i, j)
          IF (area >= threshold) THEN
             channel % mat (i,j) = 1
          ELSE
            channel % mat (i,j) = 0
          END IF
        END IF
      END DO
    END DO
  
  CASE ('ask')
     !compute slope in radians
    CALL DeriveSlope (dem, slope)
    
    DO i = 1, mask % idim
      DO j = 1, mask % jdim
        IF (mask % mat (i,j) /= mask % nodata ) THEN
           area = flowAcc % mat (i,j) * CellArea (flowAcc, i, j)
           ask = area * TAN (slope % mat (i,j)) ** scale
           IF (ask >= threshold) THEN !channel initiation found
             !set channel on all downstream cell
             iu = i
             ju = j
             DO
               channel % mat (iu,ju) = 1
               CALL DownstreamCell (iu,ju,flowDir % mat(iu,ju),id,jd)
    
               IF ( CheckOutlet (iu,ju,id,jd,flowDir) ) EXIT
               iu = id
               ju = jd                            
             END DO 
             
           END IF
          
        END IF
      END DO
    END DO 
    CALL GridDestroy (slope) 
 !other methods to be implemented
END SELECT

RETURN
END SUBROUTINE ChannelInitiation